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ABSTRACT 

We investigate whether it is possible to study perturbatively the transition 
in cosmological clustering between a single streamed flow to a multi streamed 
flow. We do this by considereing a system whose dynamics is governed by 
the Zel'dovich approximation (ZA) and calculating the evolution of the two 
point correlation function using two methods: 1. Distribution functions 2. 
Hydrodynamic equations without pressure and vorticity. The latter method 
breaks down once multistreaming occurs whereas the former does not. We 
find that the two methods give the same results to all orders in a perturbative 
expansion of the two point correlation function. We thus conclude that we 
cannot study the transition from a single stream flow to a multi-stream flow in 
a perturbative expansion. We expect this conclusion to hold even if we use the 
full gravitational dynamics (GD) instead of ZA. 

We use ZA to look at the evolution of the two point correlation function at 
large spatial separations and we find that until the onset of multi-streaming the 
evolution can be described by a diffusion process where the linear evolution at 
large scales gets modified by the rearrangement of matter on small scales. We 
compare these results with the lowest order nonlinear results from GD. We find 
that the difference is only in the numerical value of the diffusion coefficient and 
we interpret this physically. 

We also use ZA to study the induced three point correlation function. At 
the lowest order of nonlinearity we find that, as in the case of GD, the three 
point correlation does not necessarily have the hierarchical form. We also find 
that at large separations the effect of the higher order terms for the three point 
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correlatin function is very similar to that for the the two point correlation and 
in this case too the evolution can be be described in terms of a diffusion process. 

Subject headings: Galaxies: Clustering - Large Scale Structure of the Universe, 
methods: analytical. 
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1. Introduction 

The inviscid hydrodynamic equations without pressure and vorticity (referred to as 
the HD equations in the rest of this paper) are often used to describe the evolution of 
disturbances in an expanding universe filled with collisionless particles that interact only 
through Newtonian gravity. The disturbances that are usually considered are such that 
initially all the particles at any point have the same velocity i.e. it is a single streamed flow. 
Such a situation is correctly described by the HD equations. As the disturbances evolve the 
particle trajectories intersect and there are particles with different velocities at the same 
point i.e. the flow becomes multi-streamed. When this occurs the HD equations are no 
longer valid. This is because the HD equations neglect the local stress tensor associated 
with the moments of the velocity about the mean velocity at a point. 

The BBGKY hierarchy of equations obeyed by the distribution functions can be used 
instead of the HD equations. The distribution functions keep track of the position and 
velocity of the particles and these equations are valid even if multistreaming occurs. The 
question we would like to address in this paper is whether we can study the effects of 
multistreaming by using distribution functions perturbatively to follow the evolution of the 
disturbances. 

We look at the perturbative evolution of the density - density two point correlation 
function for Gaussian initial conditions in a universe with Q = 1. The perturbative 
expansion of this function using the HD equations has been studied by many authors 
(Juskewicz 1981; Vishniac 1983; Fry 1994). In a recent paper (paper II,Bharadwaj 1995) 
we have calculated the lowest order nonlinear term for the two point correlation function 
using the moments of the BBGKY hierarchy. These equations are based on the distribution 
functions and are valid even in the multi-streamed regime. The two different methods of 
calculation (HD and BBGKY) aro found to give the same result at the lowest order of 
nonlinearity, and hence, to this order, distribution functions have not been able to capture 
any effect of multi-streaming. In this paper we investigate whether by going to higher 
orders of perturbation theory we shall be able to study any effects of multi-streaming or if 
it is a limitation of perturbation theory that it cannot follow the transition from a single 
streamed flow to a multi streamed flow. 

Because of the difficulty in calculating the higher order terms in a perturbative 
treatment of gravitational dynamics (GD), we look at a simpler system where we use the 
Zel'dovich approximation (ZA, Zel'dovich 1970) to determine the motion of the particles. In 
this situation too the transition from a single streamed flow to a multi-streamed flow occurs 
and we can analyse it to see if in a perturbative calculation using distribution functions we 
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can include any effects of multi-streaming which would be missed if the HD equations were 
used instead. 

In section 2 we discuss the evolution equations. In section 3 we use distribution 
functions to calculate the evolution of the two point correlation function. In section 4 we 
do the same calculation using the HD equations and compare the result with that obtained 
in section 3. 

Bond and Couchman (1988) have studied the evolution of the two point correlation 
function using ZA and the calculation presented in section 3 is on similar lines. In a more 
recent paper Schneider and Bartlemann (1995) have studied the evolution of the power 
spectrum in ZA. For a comprehensive article on various aspects of ZA the reader is referred 
to a review by Shandarin and Zel'dovich (1989). 

In paper II we investigated the lowest order nonlinear correction (using GD) to the 
two point correlation for initial power spectra of the form P(k) oc k n at small k and an 
exponential or Gaussian cutoff at large k. We found that for < n < 3, the numerical 
results for the nonlinear correction to the two point correlation function at large x could be 
fitted by a simple formula . We also interpreted this formula in terms of a simple diffusion 
process. In section 5 of this paper we investigate the evolution of the two point correlation 
function at large separations using ZA and compare it with the results from GD. 

In section 6 we look at the evolution of the induced three point correlation function 
using ZA. This was first calculated for GD by Fry (1984) and he concluded that for power 
law initial conditions, at large separations, the three point correlation function could be 
described by the hierarchical form where it can be written in terms of products of the two 
point correlation function evaluated at the separations involved. In an earlier paper (part 
I, Bharadwaj 1994) we calculated the same quantity and found that these conclusions were 
not fully correct. We showed that the three point correlation function at some length scale, 
depends not only on the two point correlation at the same length scales but on all smaller 
scales also. As a result we found that the hierarchical form is true for only a class of initial 
conditions and that there was a class for which it did not hold. In this paper we first 
calculate the three point correlation function at the lowest order of nonlinearity for ZA and 
compare it to the results from GD. We then go on to study the effect of the higher order 
nonlinear terms at large separations. 

The calculations using ZA are valid for any value of Qq but whenever we make 
comparisons with GD it is for the specific value Q = 1. 

A similar calculation has been done by Grinstein and Wise (1987) who have studied 
the evolution of skewness of the density field averaged over a Gaussian ball. Also, Munshi 
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and Starobinsky (1994) have considered the evolution of the skewness of the density field 
for ZA and various other approximations, and Bernardeau et. al. (1993) have calculated 
the evolution of the skewness of the density field averaged over top hat filters. All of these 
calculations have been done at the lowest order of nonlinearity. 

In section 7 we present a discussion of the results obtained and the conclusions. 



2. Evolution of the distribution fuction 

The Zel'dovich approximation defines a map from the initial position af a particle to its 
position at any later instant. If x M (t) is the comoving coordinate of a particle at any time 
t, the initial instant being t , and b(t) the growing mode in the linear analysis of density 
perturbations, this map is 

x^it) = x^(t ) + b(t)Uf,. (1) 
The quantity is related to the peculiar velocity v M (t) at any instant by 

v M (t) = a(t)ic M (t) = a(t)b(t)u M (2) 

where a(t) is the scale factor. 

We consider a system of particles whose motion is governed by this mapping. This can 
be described by a distribution function f(x,u,t), where f(x,u,t)d 3 xd 3 u is the number of 
particles in the volume d 3 x around the point x and having a value of u in an interval d 3 u 
around u. 

We can see that Liouville theorem is true for the mapping defined in equation (fl|). 
Using this we can obtain the equation for the time evolution of the distribution function /, 

f(x,u,t)= f(x-b(t)u,u,t ). (3) 

We can also use equation ([!]) to obtain a differential equation for the evolution of the 
distribution function 

d d 

—f(x,u,b)+u^—f(x,u,b) = 0, (4) 

where we use the growing mode b instead of time as the evolution parameter. 

We are interested in the evolution of the statistical properties of an ensemble of such 
systems. 
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Every member of the ensemble initially has the particles uniformly distributed. Initially 
each particle can be labeled by its coordinate x^. The particles are given velocities u^(x). 
The velocity field is the gradient of a function ip(x) which for each system is a different 
realization of a Gaussian random field. It is assumed that ip is statistically homogenous and 
isotropic. The statistical properties of the ensemble are initially fully specified by the two 
point correlation of ip which is defined as (p(x) =< ip(0)ip(x) >, where the angular brackets 
< > denote ensemble averageing. 

The statistical quanitity whose evolution we shall focus on in this paper is the density 
two point correlation function £ (x, t) . This is defined by the relation 

<p> 2 (l + ax))=<p(0)p(x)> . (5) 

where p(x) is the mass density. This is just the number density of particles multiplied by 
the mass of each particle which is assumed to be the same for all the particles. 



3. The two point correlation using distribution functions 

In this section we look at the evolution of the ensemble averaged two point distribution 
functions p 2 . This is defined as 

p 2 {x\x\ u\ u 2 , t) =< f{x\u\ t)f(x 2 , u 2 , t) > . (6) 

From homogeneity and isotropy we can also say that 

p 2 (x 1 , x 2 , u 1 , u 2 , t) = p(x, u 1 , u 2 , t) (7) 

where 

X fl ~ X fJL ~ X fJL' (8) 

The density two point correlation function is related to the zeroth moment of the two point 
distribution function with respect to u. 

< p > 2 (1 + £(x,t)) = J p 2 (x,u 1 ,u 2 ,t)d 3 u 1 d 3 u 2 . (9) 
In this paper we normalize < p >= 1. 

The initial two point distribution is a Gaussian in the velocities and hence specified by 
the covariance matrix 

T;l{x) =< >(*) = / ufflp^u^u'^dWu 2 (10) 
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where a, 6 take values 1,2. The initial two point distribution function then is the Gaussian 
distribution 



p 2 (x,u ,u ,t ) 



1 



W JAT(x) 



exp 



1 



where AT(x) is the determinant of the covariance matrix. In terms of the potential <fi we 
have 

< u\u 2 v >= -d,d^(x) (12) 

and 



1. 



< K u l >= -gVV(o) V 

We use equation (j3|) to obtain the time evolution of p 2 

p 2 (x, u 1 , u 2 , t) = p 2 (x - (u 2 - u l )b{t),v}, it 2 , t )- 



(13) 



(14) 



This may also be written as 

p(x,u 1 ,u 2 ,t) = J 5 3 x - (x - (u 2 - u l )b (£)) p 2 (x ,u 1 ,u 2 ,t )d 3 x . (15) 
Using the Fourier expansion of the Dirac delta function and using equation ( |TT1) we have 

I x3 



p(x, u ,u ,t) 



2tt 



exp 



ikft (x^ 



exp 



X 



(27T)VA2V) 



: exp 



-o«( T_1 )^(^) 



d 3 kd 3 x' 



(16) 



Using this in equation @ and doing the u integrals we get 

b 2 (t) 



l+£(x,t) = [-^j /exp 



where 



ik^ (x^ x^j 
2 



exp 



k^k v F^ v {x ) 



d 3 x'd 3 k, (17) 



^(s)^ = -|V 2 0(O) V + 2d /i d li <j>(x) 



Doing the k integral we obtain the two point correlation as 
1 r 1 



l + £(x,t) 



■ exp 



2b 2 (t) 



2(,\ \ X /J, X l^) \ X l> X ") (x ) 



(2vr)5 b 3 (t) J yAF(x') 
Instead of integrating equation fllT]), if we do a Taylor expansion of 

b 2 (t) 



d 3 x . (19) 



exp 



'k^kyF^yix 
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and then do the k and the x' integrals, we obtain 



l + £(x,t) = 




Nowhere above has any assumption been made about the number of streams in the 
flow . Equation (^) obviously has the effects of multistreaming built into it. Equation ( pUD 
is what one would obtain if one did a perturbative expansion of the distribution function 
and calculated the two point correlation function. Whether by doing the perturbative 
analysis this way (i.e. using distribution functions) we are able to include the effects of 
multistreaming is what has to be checked. 



4. The two point correlation using the hydrodynamic equations 

In this section we shall work in the single stream approximation. We cosider any one 
member of the ensemble described previously. Its evoulution is described by equation (f|). 
If we take the zeroth moment of this equation with respect to u. Using the definitions 



p(x,b) = m J f(x,u,b)d u (21) 

and 

p(x, b)vn(x, b) =m J u il f(x, u, b)d 3 u (22) 
we have the continuity equation 

^p(x, b) + d„(p(z, b)v„(x, &)) = . (23) 

Next, taking the first moment of equation (§) and using equation (|23|) we have 

d 

p(x,b)[—v^(x,b) +v 1/ (x,b)d v v li (x,b)] + 

+ md u J (v u (x,b) - u u )(Vfj,(x,b) - u fl )f(x,u,b)d 3 u = 0. (24) 

In the single stream approximation the last term in the above equation is dropped, and we 
have 



d_ 

db 



v^x, b) + v v (x, b)d v Vfj,(x, b) = . (25) 
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We shall use use equations ( p3|) and to perturbatively evolve the density and velocity 
fields of the system. We then take ensemble averages and use these equations to calculate 
the two point correlation function. 

Using equation (^) we can obtain an equation for the first derivative of the two point 
correlation function 

^[< p > 2 (1 + £(x,b))] = -< dl(p(x l )v,(x l ))p(x 2 ) >-< P (x l )dl(p(x 2 )v,(x 2 )) > . (26) 
Using the normalization < p >= 1, the above equation may be written as 

J^(s,6) = -^<p(l)t;«p(2)> • (27) 

We can use equation ( |2"3"|) and (|2!S| ) to obtain equations for the higher derivatives of the two 
point correlation 

— £(*,&)) = (-ir5£fl* < P(1)«...<P(2) > . (28) 

Next we write the two point correlation function as a Taylor series in powers of the growing 
mode b 

oo in fln 

£M) = E^MW (29) 

It should be noted that this allows us to express the two point correlation function at any 
instant in terms of the derivatives of the two point correlation function at the initial instant. 
Next, using equation ( p8"D we get 

£M) = E -^r^ d ;\ d Z- d Z < p(i)«-0( 2 ) >o=o • • (so) 

n=l U - 

Then using the fact that the initial density is uniform, we can write the two point correlation 
function at any arbitrary time in terms of the initial velocities only i.e. 

oo in 

CM) = E rr(-i) n ^-^ < >*=o • • (si) 

n=l "" 

Also the initial velocity field is Gaussian and hence all the odd terms in equation (|3l|) are 
zero. We can then write this equation as 

oo i9n 

£( x b) = Y - — — d® 1 d bl . . . df" d b " < v^v b J...v^v b - > b - . . (32) 

n=l \ z ' lb )- 

For the Gaussian initial velocity field we have 

< <^<^o-< n ^ n >=Y < v%v% >< vjv h v l > ... < v^vt" > , (33) 

Ml v \ M2 ^2 Mn Vn Z— / Ml V\ fl'2 ^2 Cn fti ' \ / 

P 
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where the sum is over all possible ways of pairing the u's. 

Using this and the fact that the derivatives are symmetric in all the indices involved, 
we have for the initial velocity field 



fll U\ fJ, n lS n fll U\ fJ, n t/fi 

( ^^d b u 1 ...d a -d b u n \< v2v h v \ > ... < v?v?r > 



(34) 



This, when used in equation fl32|) , give us 



oo L2n 



>V ' / Z < ri^On Ml v l f-2 ^2 Hn 1 

n=l n - c 



< v,X >< v,X > ... < Cv > 



b=0 



(35) 



Summing the superscripts a 1; 6 1; ...a n , b n over the values 1 and 2 and using the fact that for 
the initial velocity field 



and 



we have 



< tyl >-- 



< v>l >= dffi^x) if a ± b 



-V 2 0(O) 



5^ if a = b 



(36) 
(37) 



l + £(x,t) 



oo jJ2n 



n=0 



nl 



V 2 0(O)' 



vv(oy 



(38) 



This is the same as equation (p0| ) which was obtained using distribution functions. 
Thus we see that the perturbative calculation of the two point correlation function using 
distribution functions has no effects of multistreaming and hence we reach the conclusion 
that it is not possible to perturbatively follow the transition from a single streamed flow to 
a multi streamed flow. 



5. The two point correlation at large separations. 



In this section we investigate the evolution of the two point correlation function in 
the regime where it can be studied perturbatively and we look at the behaviour at large 
separations. The initial conditions for the evolution of the cosmological correlations may be 
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expressed in terms of the potential 4>(x) or equivalently in terms of the matter two point 
correlation in the linear epoch, ^'(x, t). The two are related by the equation 

^ 1 \x,t) = b 2 {t)V*<j){x). (39) 

Usually the initial conditions are given in terms of the matter two point correlation £^\x, t) 
or its Fourier transform b 2 (t)Pi(k) which is the power spectrum. One then has to invert 
equation ( |3"9"D to obtain the potential (f)(x) and its derivatives. In doing so one has the 
freedom in choosing boundary conditions and the effect of changing the boundary condition 
is 

V 2 0(x) -> V 2 0(x) + Ci (40) 

and 

<f>{x) -> <f>{x) + ^f - + C 2 . (41) 
o 



Equation (20) for the two point correlation function is invariant uder these transformations 
and we are free to choose any boundary condition. For initial conditions where the 
integral J" °° t)xdx (or / Q °° Pi(k)dk) is finite we can choose the boundary condition 

limit a; ^ oo V 2 0(x) = 0. We then have 

roo 

< u 2 >= -V 2 0(O) = / £ {1 \x)xdx. (42) 

Jo 

In addition, if at large x the function d^d v (j){x) is monotonically decreasing and 
d,j,d u (j)(x) -C (^ V /3)V 2 0(O), we can then neglect all but one of the <9 M <9„0(a;) terms that 
appear in equation (p0f). For initial conditions where the power spectrum has the form 
P(k) oc k n at small k and if it has a cutoff at large k, the conditions discussed above are 
satisfied for n > —1. For these cases we obtain for the two point correlation function at 
large x 

5(M) .|^(=v^))" (vTvVw . (43) 

Using this we obtain the lowest order nonlinear correction to the two point correlation 
function at large x, 

e\x,t) = b ^<u 2 >V 2 e\x,t) (44) 

In paper II we have calculated the same quantity using GD and we found that for < n < 3 
at large x the results can be fitted by the formula 

e (2) (x, t) = . 1946 2 <u 2 > V 2 e (1) (x, t) (45) 



We find that the two equations are very similar and they differ only in the numerical 
coefficient. In paper II we also interpret equation fl4~5|) in terms of a simple heuristic model 
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based on a diffusion process. We consider a particular member of the enemble and look at 
the evolution of the density in volume elements located at a separation x from each other. 
We assume that the density in each volume element grows according to linear theory and 
the volume elements get rearranged randomly on small scales because of their peculiar 
velocities. Based on this model we obtained an equation identical to equation flUf) . Thus 
we see that this model gives an exact description of what happens in ZA at large scales in 
the regime when the perturbative treatment is valid. In ZA, like in our heuristic model, the 
velocity of the particles is fixed whereas in GD the particle velocity changes as evolution 
proceeds. We believe that this is responsible for the smaller diffusion coefficient for GD 
compared to ZA. 



Going back to equation Q4"5| ) and writing it in Fourier space we obtain for the power 
spectrum 

00 1 (-b 2 k 2 <u 2 > 
n! 



p(k,t) 



.n=0 



(46) 



Summing up the terms in the square brackets we have 



P(k, t) = exp 



-b 2 k 2 < u 2 > 



6 2 Pi(A;). 



(47) 



which in real space gives us 
where 



(v^2L(t))3 Jo 



cxp 



(x — x y 



^(x\t)d 3 x' 



1 



L 2 (t) = -b 2 (t) < u > . 



(4£ 



(49) 



The length scale L(t) is the r.m.s. deviation of the particles from their Lagrangian (or 
initial) positions at any time t. We see that the nonlinear evolution of the two point 
correlation function at large x corresponds to a convolution of the linear two point 
correlation with a Gaussian whose width is proportional to L(t). This is consistent with 
our interpretation of the evolution in terms of a diffusion process. 

For the case when the initial power spectrum has the form 

Pi(fc) = Ae- a2k2 k n , (50) 
using equation (^3|) at small k, we have for the nonlinear power spectrum at small k 

Pi(fc) = Ae~ {a2+L2{t))k2 k n , . (51) 
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Using equation (|50|) and (0), and using the fact that 




(52) 



we obtain for the nonlinear two point correlation function at large x 



3 + n 




2 



£i(x,t) 



1 + 




(53) 



This formula relates the nonlinear two point correlation at some separation x at a time 
t to the linear two point correlation at a smaller separation at the same time. Thus, at 
large x, for small values of the two point correlation, we have information being transferred 
out from the smaller scales to the larger scales. 

We next numerically investigate the evolution of the two point correlation function 
at large separations for the initial power spectrum P\{k) = .5e _fc2 fc. Figure 1 shows the 
function £^(x) as a function of x. This function multiplied by the square of the scale 
factor gives the linear two point correlation ^'(x,t). At large x the function has 
a negative sign and a power law behaviour x~ A . We investigate the evolution of the two 
point correlation function at the large separation x = 10. We do this using four different 
approximations which we list below: 

(a) , linear perturbation theory 

(b) . linear theory + the lowest order nonlinear correction using GD (paper II). 

(c) . the result obtained from summing the whole perturbation series for the ZA with 
the extra assumptions about the evolution at large separations made in this section i.e. 
equation flo"3"| ) 

(d) . the non-perturbative two point correlation calculated using ZA ( |19|) 

This exercise allows us to investigate two different issues. The first thing that we can 
check is how well ZA approximates GD. This can be done by comparing (b) with (c) and 
(d). In this section we have made some assumptions about the large x behaviour of the two 
point ocrrelation function and arrived at the diffusion picture for the evolution. We can 
put these assumptions to test by comparing (c) with (d). The results are shown in figure 
2. We find all the four approximations match in the early stages of the evolution. The two 
point correlation at this separation is initially negative and this value evolves according to 
linear theory where it gets multiplied by b 2 . The different approximations start to differ 
as the evolution proceeds. The first thing to note is that they start to differ much before 
£(x,t) ~ 1 when one would naively expect the perturbation series to break down. This is 
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a consequence of the non-local nature of the nonlinear terms for the two point correlation. 
As discussed in paper II, this can be understood using equation ([42]) 



which shows that the nonlinear correction depends on the linear two point correlation 
condition at all the scales and the major contribution to this integral comes from the 
small scales. The small scales become strongly nonlinear very early in the evolution and 
it is because of this that the nonlinear term starts contributing at large x even when 
£(x, t) <C 1. In all the approximations (i.e. (b),(c) and (d)) the effect of the initial deviation 
from the linear theory is to make the growth rate faster than b 2 (t). In the initial stages 
approximations (b), (c) and (d) exhibit qualitatively similar behaviour but as the evolution 
proceeds we find that (d) starts showing a behaviour completely different from (b) and (c). 
We find that the rapidly decreasing function (d) slows down its decrease and then starts 
to increase. This is quite different from the behaviour of (b) and (c) which continue to 
decrease. This difference is because of the effects of multi-streaming. In ZA the correlations 
get washed out after multi streaming occurs. Until the onset of multistreaming the diffusion 
picture (c) matches quite well with the full ZA i.e. (d). A comparison of (b),(c) and 
(d) shows that ZA qualitatively predicts the same behaviour as GD and the quantitative 
difference may be attributed to the difference in the diffusion coefficients. In the case of 
the actual gravitational dynamics (non-perturbative) we expect that the results may be 
different because there the particles will get 'stuck' in bound objects once multistreaming 
occurs (e.g. the adhesion model; Gurbatov, Saichev and Shandarin 1989) As a result of 
this the mean square displacement of the particles will be much less than in ZA or in 
perturbative GD. Although we expect this diffusion picture to hold for the actual evolution 
of the two point correlation function at large x, the perturbative treatment of GD and also 
calculations using ZA may overestimate what would be obtained in N-body simulations. 
Incidentally, the regime treated here would be difficult to study using such simulations 
since it involves the low amplitude tail of the two point correlation function which would 
be limited by the size of the box and it would require a large dynamical range. 



We use ZA to follow the evolution of the N point correlation function. It is possible to 
do this nonperturbatively by following a line of reasoning very similar to that in section 3. 
However since ZA is a good substitute for the gravitational dynamics only in the weakly 
nonlinear regime we prefer to carry out the investigation perturbatively. 




6. The 3 point correlation function. 
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We first consider the evolution of the ensemble averaged N point distribution function 
p^{x a 1 u a 1 t). This is a generalization of the ensemble averaged two point distribution 
function introduced in section 3 and the superscript a refers to the various points i.e. 1,2... 
N in phase space which are arguments of this function. Using equation (^) we obtain for 
the time evolution of this function 

p N (x a , u\ t) = p N (x a - b(t)u a , u\ t ) . (54) 

Expanding this in a perturbative series and using ai,a 2 ...a n for n indices that 
independently take values between 1 and N, and using pi, p 2 ... p n for n corresponding 
Cartesian components, we have 

oo (_U\n 

P ^x\u\t) = T,^K\Kl---K:^Z--- d Zp^ a ^ a ^)- (55) 

n=0 n - 

For both the kinds of indices the Einstein summation convention holds and all the a l s 
have to summed over the range 1 to N whenever they appear twice and the piS have to be 
summed over the three cartesian components whenever the indices are repeated. 

To calculate the N point correlation function we take velocity moments of the N point 
distribution function 



oo (_U\n 

p N (x a ,u a ,t)d 3N u = £ ±-f-d£d» < « ...nil > . 



(56) 



n=0 



All the terms where n is odd are zero and only the terms with even n contribute. We also 
have 



< u£u£ ... <:« ••• < >=< «> ... < > +permutations . (57) 

Using the fact that d^&j^ ... is symmetric in all the indices, we can add up all the 
permutations to obtain for the terms with even n 



< u a } u b } ...u an «' bn 



a U 7 >z 



2 n n\ 



(5* 



ia n b„ 



where =< u a „u v > is the covariance matrix introduced in section 3 generalized for the 
N point distribution function. 



Using this in equation (|55|) , we have 

oo (l2\n 
]3N„. _ v° 



ps(x a ,u a ,t)d^u 



oo 

On-nl Ml v\ "' Mn u v n 
n=0 Z n - 



T: 



rpa n b n 



(59) 
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In the above equation, for a fixed value of n, there will be a term with n pairs 
(cii&i), (aib2)...(a n b n ) where each index is independently summed over values 1 to N. Thus, 
for a fixed value of n, the total contribution is a sum of N 2n terms each corresponding to a 
different set of values for the position indices. In any one of these N 2n terms there can be 
two kinds of pairs 

A. if di = bi, then Tjffy = _I<5^.v 2 0(O) is a constant 

B. if di 7^ bi then Tf*%* = d?*d£i(j)(ai,bi) is a function of the separation between these two 



Any of the terms can be representd by a directed graph with N vertices and n edges. The 
pairs of the kind A correspond to an edge connecting a vertex to itself and a pair of the 
kind B corresponds to an edge connecting two different vertices (figure(3)). The integral 
J °° p^(x a ,u a ,t)d 3N u then corresponds to a sum of graphs with N vertices and the number 
of edges going from to infinity. 

The quantity / °° p^(x a , u a , t)d 3N u d 3 x 1 d 3 x 2 .. d 3 x N is the mean number of particles we 
expect to find in the volume d 3 x l at x 1 and d 3 x 2 at x 2 and ... d 3 x N around x N . This has 
contribution from the lower (i.e N-l,...l point) correlation functions also. The residue when 
the contributions from the lower correlation functions have been removed, is called the 
reduced N point correlation function. Henceforth we shall refer to the reduced N point 
correlation function as the N point correlation function. The graphs that do not connect all 
the N points correspond to functions that do not refer to all the N points and these are the 
contributions from the lower correlations. The reduced N point correlation can be calulated 
by considering only the connected graphs with N vertices. The lowest order contribution to 
the N point correlation corresponds to the connected graphs with the least number of edges. 
These graphs are the tree graphs and they have N-l edges. The other terms that contribute 
to the N point correlation can be generated by adding more edges to the tree graphs. 

We use equation ( |59|) to calculate the three point correlation function. The lowest 
order at which the three point correlation develops is n = 2 and this can be written as 



points. 




where o 1; a' 2 and a 3 are to be summed over all possible permutations of 1, 2 and 3. Equation 
( |60f) correponds to the only possible tree graph with three vertices a l: d 2 and a' 3 , and two 
edges (a'^a'2) and (a' 2 ,a 3 ) (figure(4)). 



Using 




(61) 
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we have 



C (1) (l,2,3,t) 



[i + cos 2 9 xy )e\x,t)e\y,t) 



+ cos e^—^ix, t)y^)(y t t) + ~(1 - 3 cos 2 6 xy )^\x, t) 

- l(i-3cos 2 e xy )^(x,t)^(y,t) 



(62) 



where 



and 



x — x 2 — x 3 , y — x 2 - x 



'xy 



xy 



(63) 



This explicitly exhibits the dependence of the lowest order induced three point 
correlation function on the initial two point correlation function. We see that the three 
point correlation depends on both ^'(x, t) and £^(x, t). Thus we see that the small scales 
can influence the three point correlation at large scales through the quantity £W(a;,t). The 
lowest order induced three point correlation function calculated using ZA is very similar 
to that calculated by studying gravitational dynamics perturbatively at the lowest order 
beyond the linear theory (paper I) and the difference is only in the numerical factors . 

We next calculate the higher order terms that contribute to the three point correlation 
function. These are generated by adding more edges to the tree graphs. Figures 5 and 6 
illustrates the simplest cases where we add only one edge to the tree graph. Next consider 
any of the graphs with n > 2 edges. In this graph the tree graph can be embedded in CV; 
ways. Using this in equation (R3) we have 



C(l,2,3,t) 



&2(>+2) 



oo 

On+l 
n=0 ^ 



TV. 



d ai d bl 



f)a n ab n ocij oa 2 oa 2 oa 4 



■■■■ 1 y. n u n a\ai oizola 



T, 



aibi 



(64) 



As discussed in the previous section, at large x the contributions from the terms with 
di = bi will dominate, i.e. at the lowest order, graphs of the kind shown in figure 5 . Thus, 
at large x the three point correlation function may be written as 



C(l,2,3,t) = E 



n=0 



L2n 



2 n n\ 



;-iv^(0))"(V ai ) 2 (V ai ) 2 .. (v a ") 2 d 1} (i, 2, 3, t) 



(65) 
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where the index a { indicates at which point the Laplacian acts, and it is to be summed over 
the values 1, 2 and 3. In Fourier space we have 



oo yin j 



F 3 (k\ k\ k 3 , f) = £ ^(^ 2 mn(k ai ) 2 + (k a2 ) 2 ... (*"») 2 )if '(fc 1 , k\ k 3 , t) (66) 



n=0 



2 n n\ y 3 



where F% is the Fourier transform of the three point correlation and F^ is the Fourier 
transform of the lowest order three point correlation. The terms can be summed up to 
obtain 



F 3 (k l ,k 2 ,k 3 ,t) = exp 
which gives us in real space 
C(x\a; 2 ,a; 3 ,f) = 



1 b 2 < u 2 > 



((k 1 ) 2 + (ky + (k 3 Y) 



2\2 , /;„3\2\ 



Ft\k\k 2 ,k 3 ,t) (67) 



_J / 

V^7rL(*)) 9 J 



exp 



(x a - y a f 
2L 2 {t) 



(V\y\y 2 ,y 3 ,t)d%. (68) 



Thus, at large separation, the effect of including the higher order terms for the three point 
correlation function is to convolve the lowest order induced three point correlation with 
a Gaussian of width L(t). As with the two point correlation function, this too can be 
interpreted in terms of a diffusion process. 



7. Discussion and Conclusions. 



We find that when we calculate the two point correlation function as a series in powers 
of the growing mode, we get the same answer if we do the calculation using distribution 
functions or if we do it in the single stream approximation. Since the first method is valid 
even after multi streaming occurs and the second method breaks down once multistreaming 
occurs, once multi streaming has occured we would expect to get different answers using 
the two different methods. But the two results match to all orders in the expansion 
parameter. We therefore conclude that even though these equations are valid in the multi 
streamed epoch, if we start from single streamed initial conditions we cannot perturbatively 
calculate any effect due to multistreaming e.g. vorticity, pressure. This limitation arises 
from the fact that the full two point correlation function for ZA, which includes the effects 
of multistreaming, is an exponential in p. All the derivatives of the function ^ vanish 
at b = 0. As a result, if we try to expand this function in a series in powers of b around 
b = 0, we find that coefficients of all the powers of b are zero. If one considers the power 
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spectrum instead, it is of the form e~ ak b . This function can be expressed as a power series 
in b 2 and one might that it is possible to perturbatively study the effects of multi-streaming 
by working in Fourier space instead of real space. Such a conclusion would be erronous as 
none of the terms in this expansion would have the effects of multi-streaming. It would be 
possible to study the effects of multi-streaming only if it were possible to sum the whole 
series. This point is further illustrated in an appendix where we consider a simpler example 
where a similar situation occurs. 

Shandarin and Zel'dovich (1989) present a formula for N, the mean number of streams 
at any point, in a situation where the particles are moving in one dimension under ZA. At 

-A 

small b this formula is of the form N = I + where A is a constant characterizing the 
initial conditions. If we expand this in powers of b, the coefficients for all the terms are 
zero and we find that the mean number of streams is one. This confirms that the effects of 
multi streaming cannot be studied perturbatively. Although in this analysis we used ZA, 
we expect this to hold for the full gravitational dynamics too, as derived at the lowest order 
of nonlinearity in paper II. 

In our comparison of the two point correlation function at large separations we find 
that the results obtained using ZA are quite similar to the lowest order nonlinear results 
obtained using GD and both of them can be interpreted in terms of a diffusion process 
where the rearrangement of matter on small scales affects the two point correlation at large 
scales. In ZA, for an initial power spectrum with n > — 1, the mean square displacement 
of the particles from their original positions is L 2 (t) = b 2 (t) < u 2 > and this makes its 
appearance in the formula for the nonlinear corrections to the two point correlation function 
obtained using ZA. Interpreting the results from GD in a similar fashion, for an initial 
power spectrum with n > 0, we have L 2 (t) ~ .58& 2 (t) < u 2 >. In paper II we also considered 
the case with n = and for this case we found L 2 (t) ~ \A% 2 (t) < u 2 >. The differences 
can be understood in terms of the fact that in ZA the particles move along trajectories 
calculated using linear GD, whereas when we take into account nonlinear corrections, the 
trajectories get modified by the tidal forces. In the equations for the evolution of the two 
point correlation function the tidal force acts through the three point correlation function. 
The tidal force of the third particle (in the three point correlation), will cause the other two 
particles to move towards or away from one another. This effect will be strongly dependent 
on the spatial behaviour of the three point correlation function. For the cases with n > 
the induced three point correlation has the hierarchical form at large x whereas s for the 
case with n = the induced three point correlation does not have this form. We propose 
that it is because of this that the effect of the tidal forces is different in these two cases and 
in the former case the effect of the tidal forces is to reduce the mean square diplacements 
relative to ZA whereas in the latter case it increases it. Thus indirectly, it is a diagnostic 
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of the effect of the backreaction of the three point correlation function on the pair velocity 
which in turn effects the two point correlation. 

We find that for ZA, at large x, we can sum up all terms in the perturbative series and 
the nonlinear two point correlation function is related to the linear two point correlation 
by a convolution with a Gaussian of width oc L(t). We also find that for special initial 
conditions where the power spectrum has a Gaussian cutoff at large k, the evolution at 
large x can be described by a simple scaling relation according to which the information 
propagates outward. 

We also find that this picture based on diffusion gives a good description of the 
evolution under ZA until the onset of multistreaming. Based on this we suggest that the 
evolution of of the two point correlation function in GD can also be described by a diffusion 
process until the onset of multistreaming. 

We have calculated the lowest order induced three point correlation function using ZA 
and we find that it is very similar to the result obtained using GD and the two differ only 
in the numerical factors. We also investigate the effect of the higher order nonlinear terms 
and we find that at large x we can sum the whole perturbation series. We find that the 
expression obtained after taking into account the nonlinear corrections is related to to the 
lowest order three point correlation function by a convolution with a Gaussian of width 
oc L(t). This is very similar to the evolution of the two point correlation function at large 
separations. 

It can be shown that a similar relation holds for the higher correlation functions also 
but we do not pursue this matter in this paper. 

The author would like to thank Rajaram Nityananda for his advice, encouragement, 
and for many very useful suggestions and discussions. 



Consider a Gaussian function of the variable x with standard deviation a. We are 
interested in the power series expansion of this function in a around a = 0. We can do this 
expansion by taking the Fourier transform of the Gaussian i.e. 



A. Appendix 




(Al) 
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and then doing a Taylor expansion (convergent) of e 2 fc2 °" 2 . We then get 




which gives us 

- 7 ^^ = tU 1 -a^) n S( X ). (A3) 
V27ro- ^ n\ \2 dx l ) 

Equation flX3| ) can also be derived if we take the Gaussian function and directly do a 
Taylor expansion in a i.e. without going to Fourier space. 

We see that the series expansion is entirely made up of Dirac delta functions and its 
derivatives and hence it has nonzero value only when x — 0. This should be compared with 
the original Gaussian function which has nonzero value even if x 7^ 0. We see that in this 
case the Taylor expansion fails to capture an important aspect of the original function and 
we can attribute this to the fact that we are doing the Taylor expansion of a function which 
is an exponential in 

If instead of working in real space we work in Fourier space, we find that we have to 
deal with the Taylor expansion of a function which is an exponential in a 2 instead of -j. 
There is no problem in expanding this function in a Taylor series and on might be led to 
think that the limitation of the Taylor expansion in real space can be overcome by going 
to Fourier space. But this turns out to be wrong. On comparing equations (|A2|) and ( |A3| ) 
we see that each term in the expansion in Fourier space corresponds to some derivatives 
of a Dirac delta function and hence it cannot capture any of the effects missed out if the 
analysis is done in real space. These effects can be included only if one is able to sum the 
series in Fourier space and then do the Fourier transform. 
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Fig. 1. — The initial two point correlation as a function of the separation for the power 
spectrum P(k) = .5e~ k k. 
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b(t) 



Fig. 2. — The two point correlation at a fixed separation x = 10 as a function of the growing 
mode b(t) for (a) linear theory, (b) linear theory + lowest order nonlinear correction using 
GD (c) nonlinear evolution using ZA and the assumptions made in section 5 about the large 
x behaviour, and (d) nonperturbative ZA 
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Fig. 3. — This shows the two possible kinds of edges A. connects a vertex to itself B. connects 
two different vertices. 
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Fig. 4. — This shows the tree graph corresponding to the lowest order induced three point 
correlation function. 



Fig. 5. — This shows some of the graphs corresponding to the contribution to the three point 
correlation function at one order beyond the lowest. These graphs are all obtained by adding 
edges to the tree graph. These graphs show those cases where the extra edge connects a 
vertex to itself. 
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Fig. 6. — This shows some of the graphs corresponding to the contribution to the three point 
correlation function at one order beyond the lowest. These graphs are all obtained by adding 
edges to the tree graph. Thee graphs show those cases where the extra edge connects two 
different vertices. 



